Notes on zoomed in maps and related topics

You can download the materials for this notebook from this link.

First load the usual libraries

library(sf)

## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.1.0; sf_use_s2() is TRUE

library(tmap)

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

library(tmaptools)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Next, we’ll load a couple of datasets. These are from the week 7 lab materials.

abb <- st_read("abb-2193.gpkg")

## Reading layer `abb-2193' from data source 
##   `/home/osullid3/Documents/teaching/Geog315/slides/zoomed-in-maps/abb-2193.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2125 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1742685 ymin: 5420357 xmax: 1785682 ymax: 5492218
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000

welly <- st_read("wellington-base-data.gpkg")

## Reading layer `wellington-base-data' from data source 
##   `/home/osullid3/Documents/teaching/Geog315/slides/zoomed-in-maps/wellington-base-data.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 122 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1735093 ymin: 5410973 xmax: 1774833 ymax: 5444626
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000

These layers are in the same projection, so we don’t need to worry about that, but as is usual, it is a good idea to transform all layers to a matching projecting suitable for the project. Use st_transform() to do this.

Web maps

The easiest way to handle variations in scale and extent of layers is a web map. For this, use tmap_mode("view")

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(welly) +
  tm_polygons(col = "pop", alpha = 0.5) + 
  tm_shape(abb) +
  tm_dots()

The different extents of these two datasets might be a problem. If we are OK with filtering the point data to match the polygons data, then just use st_filter():

welly_abb <- abb %>%
  st_filter(welly)
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(welly) +
  tm_polygons(col = "pop", alpha = 0.5) + 
  tm_shape(welly_abb) +
  tm_dots()

With a web map you can use the export options in RStudio to make snapshots of zoomed in portions of your maps for inclusion in a report.

When a web map is no good

Sometimes a web map is not a viable option. We have a number of static map options available. Before looking at them, I will switch back to static mapping:

tmap_mode("plot")

## tmap mode set to plotting

Bounding box

We can limit the extent of a static map with a bbox option in the first layer. By default, tmap uses the extent of the first layer to set the bounding box. So

tm_shape(welly) + 
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

makes a map that only extends as far as that first layer. If we want it to extend to the full extent of the abb data then we can do

tm_shape(welly, bbox = abb) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

We’ll worry about the missing ‘base’ layer later in this document!

We can also use a bounding box to zoom in on an area. For example, if we make a bounding box that is from only the Lower Hutt City TA:

lower_hutt <- welly %>%
  filter(ta_name == "Lower Hutt City")
tm_shape(welly, bbox = lower_hutt) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

Again, there is missing basemap, but we’ll get to that in a bit.

Using a bounding box for a detailed local map

If we want to make a specific bounding box, then we need to figure out the relevant coordinates. For this purpose, it is useful to make a temporary map with a grid.

tm_shape(welly) +
  tm_polygons() +
  tm_grid()

Now say we wanted a map showing only the central city. This has an extent from 1,745,000 to 1,750,000 EW and 5,425,000 to 5,430,000 NS. We can use this information to make a bounding box using the bb function from tmaptools:

cbd <- bb(xlim = c(1745000, 1750000), ylim = c(5425000, 5430000))

and now we can do

tm_shape(welly, bbox = cbd) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

There are other options in the bb function, that allow you to specify it by using cx and cy (x and y centre coordinates) and width and height. The units are always the map projection units, which here are metres, so for example another way to make the above map would be

cbd <- bb(cx = 1747500, cy = 5427500, width = 5000, height = 5000)
tm_shape(welly, bbox = cbd) +
  tm_polygons() +
  tm_shape(abb) + 
  tm_dots()

This approach is a bit easier to recalculate if you want to zoom in a bit closer or zoom out a bit, by just changing the width and height of the bounding box. There is also an ext option in the bb function to increase or shrink the bounding box relative to the requested width and height.

But what I really want is a basemap

Often what we really want a web map for is the basemap.

We can get basemaps in static maps using the rget_tiles() function in the package maptiles and including it as a tm_rgb layer in the map. Here’s how we use it. First we grab the map tiles:

library(maptiles) # you may have to install this
basemap <- get_tiles(abb)

Then we make a map

tm_shape(basemap) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

If the basemap is a bit blurry, then you can zoom in a bit. You have to be careful about this, as if you increase the zoom too much it will take ages to download (and you won’t be able to read the detail anyway). By experiment zoom levels in the range 9 to 11 are about right for the Wellington region, and get_tiles is smart enough to get it about right most of the time, without you specifying a zoom level. The only other thing to watch out for is that maptiles delivers whole web map tiles, so the cropping might be a bit strange. You can fix this by setting the crop = TRUE parameter in the get_tiles() function, or just using bbox in tm_shape() as we have seen already:

# basemap <- get_tiles(abb, crop = TRUE) # this would crop the tiles to the data
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

If the colours are distracting, you can wash them out with the saturation option in tm_rgb

tm_shape(basemap, bbox = abb) + 
  tm_rgb(saturation = 0.35) + 
  tm_shape(abb) + 
  tm_dots()

Other basemaps are available than. You can see the options in the get_tiles() function help. For example, this is a nice one (if a bit distracting):

basemap <- get_tiles(abb, provider = "Stamen.TerrainBackground")
tm_shape(basemap, bbox = abb) + 
  tm_rgb() + 
  tm_shape(abb) + 
  tm_dots()

Remember that if you are overlaying polygons on a basemap you will probably want to set alpha (the opacity) to less than 1. So here is a final example.

basemap <- get_tiles(welly)
tm_shape(basemap, bbox = welly) +
  tm_rgb(saturation = 0.3) +
  tm_shape(welly) + 
  tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
  tm_shape(abb) + 
  tm_dots() +
  tm_layout(legend.outside = TRUE)

No really… here’s a final little extra. The cbd bounding box might be something we want to use in a map like this. Unfortunately a bb object is not able to be used directly to specify the limits of a web map. We have to convert it to a simple features object with a known projection. Here’s how we can do that.

cbd_sf <- cbd %>% 
  st_as_sfc() %>% 
  st_set_crs(2193) # this is the known projection in this case
basemap <- get_tiles(cbd_sf, zoom = 13) # I experimented to get best zoom
tm_shape(basemap, bbox = cbd) +
  tm_rgb(saturation = 0.3) +
  tm_shape(welly) + 
  tm_polygons(col = "pop", palette = "viridis", alpha = 0.5) +
  tm_shape(abb) + 
  tm_dots() +
  tm_layout(legend.outside = TRUE)